home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 12 / Mac Magazin and MacEasy Magazine CD - Issue 12.iso / Sharewarebibliothek / Anwendungen / Wissenschaft & Technik / Yorick / yorick11-ppc folder / include / fitrat.i < prev    next >
Text File  |  1995-04-03  |  4KB  |  150 lines
  1. /*
  2.    FITRAT.I
  3.    Polynomial and rational function interpolators.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func fitpol(y, x, xp, keep=)
  11. /* DOCUMENT yp= fitpol(y, x, xp)
  12.        -or- yp= fitpol(y, x, xp, keep=1)
  13.      is an interpolation routine similar to interp, except that fitpol
  14.      returns the polynomial of degree numberof(X)-1 which passes through
  15.      the given points (X,Y), evaluated at the requested points XP.
  16.  
  17.      The X must either increase or decrease monotonically.
  18.  
  19.      If the KEEP keyword is present and non-zero, the external variable
  20.      yperr will contain a list of error estimates for the returned values
  21.      yp on exit.
  22.  
  23.      The algorithm is taken from Numerical Recipes (Press, et. al.,
  24.      Cambridge University Press, 1988); it is called Neville's algorithm.
  25.      The rational function interpolator fitrat is better for "typical"
  26.      functions.  The Yorick implementaion requires numberof(X)*numberof(XP)
  27.      temporary arrays, so the X and Y arrays should be reasonably small.
  28.  
  29.    SEE ALSO: fitrat, interp
  30.  */
  31. {
  32.   /* find lower and upper indices in x of the interval containing xp */
  33.   nx= numberof(x);
  34.   u= digitize(xp, x);
  35.   l= max(u-1, 1);
  36.   u= min(u, nx);
  37.  
  38.   /* initialize yp to the y(x) closest to xp */
  39.   mask= (abs(x(l)-xp)<=abs(x(u)-xp));
  40.   ns= l*mask + u*(!mask);
  41.   yp= y(ns);
  42.   ns--;
  43.  
  44.   c= d= y= array(double(y), dimsof(xp));
  45.   dx= x - xp(-,);
  46.   if (dimsof(ns)(1)) {
  47.     is= ns;
  48.     is(*)= nx*indgen(0:numberof(xp)-1);
  49.   } else {
  50.     is= 0;
  51.   }
  52.   for (m=1 ; m<nx ; m++) {
  53.     ho= dx(1:nx-m,..);
  54.     hp= dx(1+m:nx,..);
  55.     den= (c(2:nx-m+1,..)-d(1:nx-m,..))/(ho-hp);
  56.     d(1:nx-m,..)= hp*den;
  57.     c(1:nx-m,..)= ho*den;
  58.     mask= (2*ns>=(nx-m));
  59.     y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask); /* this is error estimate */
  60.     yp+= y;
  61.     ns-= mask;
  62.   }
  63.  
  64.   if (keep) {
  65.     extern yperr;
  66.     yperr= y;
  67.   }
  68.  
  69.   return yp;
  70. }
  71.  
  72. func fitrat(y, x, xp, keep=)
  73. /* DOCUMENT yp= fitrat(y, x, xp)
  74.        -or- yp= fitrat(y, x, xp, keep=1)
  75.      is an interpolation routine similar to interp, except that fitpol
  76.      returns the diagonal rational function of degree numberof(X)-1 which
  77.      passes through the given points (X,Y), evaluated at the requested
  78.      points XP.  (The numerator and denominator polynomials have equal
  79.      degree, or the denominator has one larger degree.)
  80.  
  81.      The X must either increase or decrease monotonically.
  82.  
  83.      If the KEEP keyword is present and non-zero, the external variable
  84.      yperr will contain a list of error estimates for the returned values
  85.      yp on exit.
  86.  
  87.      The algorithm is taken from Numerical Recipes (Press, et. al.,
  88.      Cambridge University Press, 1988); it is called the Bulirsch-Stoer
  89.      algorithm.  The Yorick implementaion requires numberof(X)*numberof(XP)
  90.      temporary arrays, so the X and Y arrays should be reasonably small.
  91.  
  92.    SEE ALSO: fitpol, interp
  93.  */
  94. {
  95.   /* find lower and upper indices in x of the interval containing xp */
  96.   nx= numberof(x);
  97.   u= digitize(xp, x);
  98.   l= max(u-1, 1);
  99.   u= min(u, nx);
  100.  
  101.   /* initialize yp to the y(x) closest to xp */
  102.   mask= (abs(x(l)-xp)<=abs(x(u)-xp));
  103.   ns= l*mask + u*(!mask);
  104.   yp= y(ns);
  105.  
  106.   exact= (x(ns)==xp);
  107.   list= where(exact);
  108.   if (numberof(list)) {
  109.     yexact= yp(list);
  110.     yerr= 0.0*yexact;
  111.   }
  112.  
  113.   list= where(!exact);
  114.   if (numberof(list)) {
  115.     yp= yp(list);
  116.     xp= xp(list);
  117.     ns= ns(list);
  118.  
  119.     c= d= y= array(double(y), numberof(xp));
  120.     d+= fitrat_tiny;
  121.     dx= x - xp(-,);
  122.     ns--;
  123.     is= nx*indgen(0:numberof(xp)-1);
  124.     for (m=1 ; m<nx ; m++) {
  125.       di= d(1:nx-m,);
  126.       ci= c(2:nx-m+1,);
  127.       t= dx(1:nx-m,)*di/dx(1+m:nx,);
  128.       dd= (ci-di)/(t-ci);
  129.       d(1:nx-m,)= ci*dd;
  130.       c(1:nx-m,)= t*dd;
  131.       mask= (2*ns>=(nx-m));
  132.       y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask);
  133.       yp+= y;
  134.       ns-= mask;
  135.     }
  136.  
  137.   } else {
  138.     yp= [];
  139.   }
  140.  
  141.   if (keep) {
  142.     extern yperr;
  143.     yperr= merge(yerr, y, exact);
  144.   }
  145.  
  146.   return merge(yexact, yp, exact);
  147. }
  148.  
  149. fitrat_tiny= 1.e-30;
  150.